Systematic time-scale-bridging molecular dynamics applied to flowing polymer melts 
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We present a novel thermodynamically guided, low-noise, time-scale bridging, and pertinently efficient strat- 
egy for the dynamic simulation of microscopic models for complex fluids. The systematic coarse-graining 
method is exemplified for low-molecular polymeric systems subjected to homogeneous flow fields. We use 
established concepts of nonequilibrium thermodynamics and an alternating Monte-Carlo-molecular dynamics 
iteration scheme in order to obtain the model equations for the slow variables. For chosen flow situations of 
interest, the established model predicts structural as well as material functions beyond the regime of linear re- 
sponse. As a by-product, we present the first steady state equibiaxial simulation results for polymer melts. The 
method is simple to implement and allows for the calculation of time-dependent behavior through quantities 
readily available from the nonequilibrium steady states. 
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I. INTRODUCTION 

Systematic bridging the time- and length-scale gap between 
microscopic and macroscopic levels of description is "of the 
greatest importance in theoretical science" [1]. In many cases, 
this challenging task can neither be solved purely analytically 
nor by brute force computer simulations alone. This is true 
in particular for soft condensed matter like colloids, poly- 
mers, liquid crystals, with their internal structure leading to 
additional length and time scales, intermediate between mi- 
croscopic and macroscopic scales [2]. 

In recent years, effective interactions for coarse-grained 
models of soft matter systems have been derived from inver- 
sion procedures that are designed to reproduce chosen pair 
correlation functions [3-7]. While the inversion procedures 
often reproduce the static structure rather accurately, their 
naive extension to dynamical phenomena clearly failed [4]. 
This deficiency calls for a systematic approach that bridges 
simultaneously the time- and length-scale gap between two 
levels. For comparatively simple two-dimensional crystalline 
solids, a simultaneous space/time coarse-graining procedure 
was proposed recently in [8] based on renormalization group 
techniques. There, temporal coarse graining is coupled via 
the dynamical critical exponent to the degree of spatial coarse 
graining. This approach is unfortunately not applicable to 
the dynamics of complex fluids, since their internal struc- 
tures break the scale invariance - an essential prerequisite for 
renormalization group methods - and lead to the emergence of 
slow, non-hydrodynamic modes. The latter are typically de- 
scribed on an intermediate, mesoscopic level by a set of "col- 
lective" or "structural" variables II(z) which in turn deter- 
mine the macroscopic properties of complex fluids [2]. Since 
many microstates z are compatible with given values of II, 
the mesoscopic level is necessarily stochastic in nature. Thus, 
the emergence of entropy and irreversibility from reversible 
dynamics is the hallmark of coarse graining. Several coarse- 
graining approaches, in particular for solutions and suspen- 
sions, have been suggested where the starting level is already 
dissipative (see e.g. [7, 9] and references therein). In the 
context of polymer melts, promising work on coarse-graining 
polymer chains starting from Hamiltonian dynamics has been 
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FIG. 1: (color online). Components of gyration tensor x (left) and 
Lagrange multiplier A (right panel) vs. shear rate for a FENE poly- 
mer melt (iV = 20). Lagrange multipliers self-consistently enter 
the anisotropy and stretching of polymer chains. Comparison with 
standard NEMD reference results (left panel) show that the gener- 
alized canonical distribution (3) provides a good description of the 
nonequilibrium stationary state in shear flow. We use Lennard- Jones 
units throughout this paper. 



done e.g. in [10, 11]. 

In this paper, we propose and explore a systematic, thermo- 
dynamically guided method which establishes the mesoscopic 
model from the underlying microscopic level. The proposed 
method is general enough to be applied to various soft matter 
systems and valid in equilibrium as well as nonequilibrium 
situations. Its strategy relies on the balance of reversible and 
irreversible contributions to the dynamics and explicitly ac- 
counts for the entropy generated in the coarse-graining step 
[12]. We use an alternating Monte-Carlo (MC) and molecu- 
lar dynamics (MD) simulation scheme in order to iteratively 
determine static and dynamic "building blocks" [13] of the 
mesoscopic model self-consistently. 



II. ORIGINAL AND COARSE-GRAINED MODEL SYSTEM 

The novel algorithm is applicable to a wide range of soft 
matter systems. In order to illustrate the basic idea and its 
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worked out counterpart, let us consider a particular liquid, a 
classical monodisperse bulk model polymer melt. The system 
consists of N c h anharmonic multibead-spring (FENE) chains 
made of N purely repulsive Lennard-Jones beads each [14- 
16]; A^b = iV c h./V particle positions and momenta are denoted 
as {r_j} and {pj}. The interaction energy between particle i 
and j is Uij 
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(1) 

The distance between particles i and j is 
denoted by ry, a the bead diameter and e the Lennard-Jones 
interaction energy. Chain connectivity is ensured by FENE 
springs that act between adjacent neighbors along the chain, 
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and t/ EENE = for all other particle pairs. All model param- 
eters and thermodynamic state point are adopted from [14]: 
temperature T = e/k-Q, density n = 0.84 er -3 , finite exten- 
sibility of the springs ro = 1.5 a, and the strength of chain 
potential 6fene = 67.5 e is large enough in order to prevent 
chain crossings. In the following, we use reduced Lennard- 
Jones units throughout [17]. 

The simple FENE model system is very useful to describe 
the general dynamical behavior of polymer melts [14-16, 18]. 
This system serves as our starting point, providing the mi- 
croscopic ("atomistic") level of description without any ir- 
reversibility built in. Under the assumption that the collec- 
tive variables II capture all relevant physical processes on the 
time scale of interest, the nonequilibrium state of the system 
is characterized by the generalized canonical ensemble, 



p(z) = / eq (z) e 



-A:n(z)-A 



(3) 



with phase space coordinates z = {rj,pj} and the classi- 
cal /eq(z) oc exp{— H {2.) / k-^T} with H denoting the mi- 
croscopic Hamiltonian [1, 13, 19]. The Lagrange multipliers 
A(x) (cf. Fig. 1) are determined by the values of the slow 
variables, x = (II (z)), where the average is performed with 
(3), and Ao a normalization constant. As structural variable, 
we here choose x to be the mean tensor of gyration, 



n(z) = 
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where i a = (a — 1)N and r° = N^ 1 Y^=4 +i Vi * s *- ne cen " 
ter of mass of chain a. This choice of slow variables is ap- 
propriate for low-molecular, unentangled polymeric systems, 
where x indeed varies slowly compared with fast relaxation 
processes such as fluctuation of bond lengths and angles, in- 
termolecular distances, or higher normal modes [2, 13]. For 
a more detailed justification of our choice of x see Appendix 
A. We can neglect the macroscopic hydrodynamic velocity 
field in (3) since it equilibrates extremely rapidly on length 



scales of individual polymers [20]. The same situation is en- 
countered in other complex fluids as long as the large relax- 
ation time scales of the collective variables are generated on 
relatively short length scales. For a more complete treatment 
including the hydrodynamic fields see Ref. [21]. 

The time evolution for the slow variables x can in general 
be written as [ 1 3] 



x = x rcv + M : — 
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— = fc B A, 
ox 



(5) 



where x rev denotes the reversible contribution in terms of a 
Poisson bracket. Here, we have employed the expression for 
the macroscopic entropy S(x) = —fee (In p) corresponding to 
the ensemble (3). Entropy gradients drive the irreversible con- 
tribution to (5). Equation (5) is justified e.g. from projection 
operator derivation [13, 22], which shows that the symmet- 
ric friction matrix M(x) can be obtained from a Green-Kubo 
type formula 



M = (M(z(t))), M 
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■A Ts II(z)A Ts II(z), (6) 



where A T5 II denotes fast fluctuations of II on a time scale r s 
that separates the evolution of the slow variables x from the 
rapid dynamics of the remaining degrees of freedom. 

The reversible part of motion is obtained analytically by 
considering the transformation behavior of II, cf. [13] for 
worked out examples. Specifically, when x is a conformation 
tensor such as the tensor of gyration, and considering a macro- 
scopic flow field v(r) = k ■ r, hence k = (Vv) T , one finds 
the so-called upper-convected behavior [21], x rev (x, k) = 
x • k t + k ■ x. The remaining building blocks A and M 
needed to complete the coarse-grained model (5), we obtain 
self-consistently through a hybrid iteration scheme, as de- 
scribed next. 



III. SYSTEMATIC TIME-SCALE BRIDGING METHOD 

In general, the space of admissible values for the slow vari- 
able x is too large for a full parameterization of A(x) from 
direct numerical integration. We choose to parameterize A 
and M along one-dimensional paths x(7), where 7 denotes 
the value of the external control parameter, i.e. the flow rate 
for chosen velocity gradients ^(7) in our case. Note, that this 
procedure is analogous to the experimental determination of 
rheological properties in viscometric flows [2]. While errors 
in determining A can in principle violate the thermodynamic 
integrability condition for S*(x), this problem is avoided when 
working with one-dimensional paths which do not cross. In 
order to calculate A(x) for relevant x (here, relevant for given 
flow gradient k), we investigate nonequilibrium steady states, 
for which the left hand side of (5) vanishes. The systematic 
time-scale bridging method we propose is summarized in Tab. 
I. 

The updated Lagrange multipliers obtained in step (v) can 
potentially be used to re-enter the procedure at (i), and follow 
steps (ii)-(v) until A has converged. The whole procedure 
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step description 

(i) choose initial values for the Lagrange multipliers A 

(ii) generate n independent configurations distributed 
according to the generalized canonical ensemble (3) 

(iii) solve Hamilton's unconstrained equations of motion 
for all n systems during a short time interval r s 

(iv) calculate the friction matrix M from Eq. (6) and 
x directly from the n trajectories produced in (iii) 

(v) calculate an updated value for A by solving (5) 
for A with x = (in terms of M, x, and k 
the latter two quantities are "hidden" in x rC v) 



TABLE I: Summary of proposed time-scale bridging method. 



(i)-(v) is then repeated for other choices of the control param- 
eter 7 in order to establish the model (5) for different external 
fields. 

Notice, that the strategy does not require the implemen- 
tation of flow-specific boundary conditions such as Lees- 
Edwards (shear) [14] or Kraynik-Reinelt (planar elongational 
flow) [18] which is a particularly useful feature as it allows us 
to study arbitrary flow situations within exactly the same ap- 
proach. In the same spirit, and in order to not potentially fal- 
sify results for the friction matrix, the algorithm also does not 
involve any constraints such as thermo- or barostats. These 
advantages are build in our approach since the macroscopic 
variables do not change significantly on the short time scale 
r s of the MD simulations in (iii). 

We now specify how to implement the steps (i)-(v) effi- 
ciently, and how to self-consistently determine the range of 
validity of the underlying assumption (3). We choose the 
control parameter 7 logarithmically equidistant, log(7) 6 
[a, a + Aa, a + 2Aa, .., b}. Before we start the procedure, 
we initialize A = and log 7 = a. 

The loop starts at (i) with the current value of A. For (ii) 
the same A is used in a MC scheme to generate microscopic 
configurations distributed according to (3). We have gener- 
ated n realizations (typically, n = 500) by slightly modi- 
fying the procedure of [23]: For each realization, we gen- 
erate N{. h > A^ch (infinitely thin) independent single FENE 
polymer chains, each distributed according to exp(— A : II*), 
where II* = Yl/N c ^ is the tensor of gyration of the sin- 
gle chain. Next, the diameter of chains is successively in- 
creased, and overlapping chains selectively removed. With 
this method, we generate a polymer melt at the desired den- 
sity, where the anisotropy generated by A remains preserved. 
Subsequently, Maxwellian distributed velocities are assigned, 
in agreement with (3). For (iii) one chooses a symplectic in- 
tegrator (we have used a velocity-Verlet algorithm) to per- 
form microcanonical equilibrium MD based on the micro- 
scopic Hamiltonian H(z). We calculate and store trajectories 
z(i) during a short time interval, t G [0,r s ], which is small 
enough to not significantly alter x during the course of the 
MD. For polymeric systems, the gyration tensor will relax to- 
wards equilibrium on a time scale r which is known to be huge 
compared with the Lennard-Jones time unit, r = 0.39 (1 + 



N/78)N 2 from [14] for melts under study, i.e., r w 200 for 
N = 20. As we carefully investigated, results are (as they 
should for proper choice of r s ) insensitive on r s in the regime 
t/t s e [5,50]. See also Fig. 2. We use t s = t/30 <C t, 
and N G {10,20,30} for results to be presented. Notice, 
that the MD simulation time is thus very short compared to 
conventional nonequilibrium MD (NEMD) at (the problem- 
atic) low field strengths (flow rates), where simulation times 
large compared with the inverse rate (7 -1 ) are required, (iv) 
With the n sets of phase space trajectories z(t) at hand, one 
inserts them into the definition of the slow variable n(z(t)), 
and then evaluates the friction term M (in our case a 4 x 4 ma- 
trix) from (6), with A Ts LT(z) = LT*(z(t s )) - II*(z(0)). The 
average in (6) is evaluated as an arithmetic mean over the n in- 
dependent trajectories, e.g., M = (1/n) Y^i -M-u), where we 
denote the partial contribution from trajectory i £ {l,..,n} 
by a bracketed subscript. The number of samples n has to 
be chosen large enough to calculate M sufficiently accurate. 
In our case, several components of M should vanish by sym- 
metry consideration, and one can choose n as large as to en- 
sure these components vanish within statistical uncertainty. 
Notice further, that M possesses basic symmetries such as 
Map^v = = Mp ajJjV for arbitrary choices of indices 

because II is symmetric, (v) Repeating the procedure (i)- 
(v)-(i)-.. for each 7 until convergence can be replaced by 
an efficient reweighting scheme. This scheme relies on the 
smallness of the change of increment Aa, which comes to- 
gether with moderate changes of the distribution function p. 
To this end we use Broyden's method with standard settings 
[24] which does not require the Jacobian matrix, to solve the 
nonlinear system 

" e -,5A:n (l) 
= [C< + k B M (i y. SA] Wi , Wi = -tA-.n,* ^ 

i=i ^3 e 

for (matrix) <5A, with mismatch C, = x rev (n(j), k) + 
fceMji): A, cf. Eqs. (3), (5). For example, in a shear flow, 
(7) stands for six equations and six unknowns. With the so- 
lution 5 A of (7) at hand, we directly calculate the reweighted 
slow variables and friction matrix, x = J2iWiT^-(i), M = 
J2i WiAt(i), as well as updated Lagrange multipliers, A — > 
A + SA. A justification of the reweighting scheme is given 
in Appendix B. Finally, we increase the control parameter 
log 7 — > log 7 + Aa, and start over with step (i) of the proce- 
dure, until we have swept through the control parameter space. 

By then, we have recorded consistent sets x, M, as well as 
A for the whole range of parameters 7. That is, we have ob- 
tained A(x) and M(x) and therefore established the coarse- 
grained model (5) for particular parameterized path x(7). By 
choosing the control parameters appropriately, our approach 
uses paths to explore those regions in state and parameter 
space that correspond to driven nonequilibrium situations of 
interest. For the system under study, the quantity x(7) itself 
is experimentally accessible by means of small angle neu- 
tron scattering [14]. Other particularly interesting material 
functions are flow curves, i.e., stress tensor er as function of 
the control parameter 7. The macroscopic expression for the 
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polymer contribution to the stress tensor 

a = -2n p fc B Tx -A, (8) 

where n p is the polymer concentration, follows from both, 
nonequilibrium thermodynamics [13], and by evaluating the 
microscopic expression for the stress tensor in the ensemble 
(3), see Appendix C. A more detailed discussion of the stress 
tensor within this context is given in [21]. 

Before presenting results obtained with the proposed 
method, we briefly comment on the time- and length-scales 
involved, already alluded to in the introduction. The origi- 
nal, microscopic model has as characteristic length scale the 
bead diameter a and reference time tlj = [ma 2 /e] 1 / 2 , where 
m is the mass and e the characteristic Lennard-Jones interac- 
tion energy. On the coarse-grained level, the characteristic 
length scale is the radius of gyration, R g « aN 1 / 2 . The 
corresponding time scale estimated from the Rouse model 
[25] is t r = ((Na) 2 /[3TT 2 k B T], where ( is the bead fric- 
tion coefficient. Therefore, the bridging of length scale 
Rg/u = N 1 / 2 is associated with a bridging of time scales 
tr/t lj = cN 2 , wherec = 5/(16TT^ 2 )[C/Co}^/k B T} 1/2 with 
Co = 37r/ (lQa^itmk-BT] 1 / 2 the friction coefficient of a hard 
sphere gas. 

IV. RESULTS 

Figure 2 shows different components of the M matrix (6) as 
a function of the separating time r s . As mentioned before, the 
results for M are to a good approximation independent of the 
precise value of t s in a broad range r s G [5, 50] which is sig- 
nificant smaller than the polymer relaxation time r (r « 200 
for N = 20). Furthermore, the comparison in Fig. 2 shows 
that the simplified formula (6) approximates the more accu- 
rate integral formula [13] M = ^ J^"dt{tl(t)tl(0)} quite 
well. 

Having established the thermodynamic building blocks 
A(x) and M(x), we can use the evolution equations (5) to 
study time-dependent flows. We have calculated transient dy- 
namics in startup of steady shear flow, or storage and loss 
moduli G" and G" as function of frequency u) upon using an 
oscillating control parameter 7 oc sinojt (graphs not shown). 
We note that, due to our choice of the parameterization x(-y), 
the transient dynamics x(i) is readily calculated as long as we 
do not leave the known subspace {x(7)}. Otherwise, interpo- 
lation and extrapolation methods are needed for parameteriz- 
ing the missing regions in x-space. 

There are several options to test the range of validity of 
the coarse-grained model. As an internal consistency check, 
we recommend comparing the macroscopic expression for the 
stress tensor Eq. (8) with the standard microscopic (virial) ex- 
pression, Eq. (CI). Both are available during the course of 
the simulation. We have verified that the two expressions for 
<7 agree with each other for the range of flow rates consid- 
ered. Under strong flow conditions and beyond the scope of 
the present study, higher order modes and kinetic contribu- 
tions to the stress tensor tend to become increasingly impor- 
tant and need to be included suitably in x, cf. [19]. 
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FIG. 2: Different components of the friction matrix M as a function 
of the separating time t 3 obtained in step (iii) of our procedure. Solid 
and open symbols correspond to the integral formula mentioned in 
the text and Eq. (6), respectively. Results correspond to a chain 
length of N = 30 and a planar shear flow with dimensionless shear 
rate 7 = 0.00036. 
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FIG. 3: (color online). Flow alignment angle \ calculated from gy- 
ration tensor x (open circles) and A (filled squares) as a function of 
the logarithm of the shear rate 7 for planar shear flow (chain length 
N = 30). Diamonds correspond to NEMD reference results taken 
from [14] obtained under the same conditions. 



We apply the proposed method to the FENE polymer melt 
described above, subjected to various flows (results for mixed 
and elliptical elongational flow not shown). For the case of 
simple shear, Fig. 1 shows the shear rate dependence of the 
chosen slow variable x (tensor of gyration in our case) and 
the corresponding Lagrange multiplier A. Very good agree- 
ment of x with NEMD reference results is obtained. As a 
further consistency check, we have verified, that the basic 
identity [x\x — 2:22)2:12 = (An — A22)Aj~ 2 1 , derived from 
Eq. (5) [26] using our choice for x and k, holds within er- 
ror margins. This quantity is related to the flow alignment 
angle x by (in — 2:22)2;^ = 2cot(2%). Therefore, we 
show in Fig. 3 the alignment angle \ calculated from x as 
well as from A. The very good agreement between those val- 
ues shows the intrinsic consistency of our results. Further- 
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FIG. 4: (color online), (a) We follow the notation employed 
in [27] where the transposed flow gradient is written as Vv = 
e [eiei + me2e2 — (1 + ra)e3e3] + 7eie2, with shear rate 7, elon- 
gation rate e, special cases m — —0.5 (simple), (planar), +0.5 
(elliptical), and 1 (equibiaxial elongation) when 7 = 0, and sim- 
ple shear, when e = 0. Besides shear viscosity 77, the graph 
shows the properly (cf. text part and [15, 27]) scaled viscosities 
Hi ee 77i/[2(2 + m)] and H 2 = r? 2 /[2(l + 2m)] vs. flow rate for 
N = 20, where 771 ee (<Tn - a: j:j )/e and 772 ee (<T22 - cr33)/e- (b) 
Maximum component of the gyration tensor xn for the same types 
of flow, vs. flow rate (N = 20). 
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FIG. 5: (color online). Polymer contribution to non-newtonian shear 
viscosity vs. shear rate for various molecular weights. Exemplar- 
ily, reference results obtained via direct NEMD simulation [14] are 
shown for N = 30. The inset shows zero-rate shear viscosity 770 and 
first viscometric function ^1.0 vs. chain length N, both coinciding 
with data from extensive NEMD [14]. 



more, our results are also in good agreement with standard 
NEMD simulation also displayed in Fig. 3 for planar shear 
flow with N = 30 [14]. Figure 4a shows shear and exten- 
sional viscosities for different flow conditions. Our results 
confirm expectations from a retarded motion expansion anal- 
ysis for a comparable system, studied via extensive NEMD in 
[15]. In particular, Fig. 4a shows that the scaled viscosities all 
superimpose for vanishing flow rates, in agreement with pre- 
dictions from linear viscoelasticity theory. Also in agreement 
with previous results, the viscosity in simple elongation ex- 
hibits a maximum around a dimensionless rate of order unity, 



while in planar and equibiaxial elongation as well as in pla- 
nar shear flow the viscosity decreases monotonically with flow 
rate [15, 27]. The corresponding x 11 -components of the gyra- 
tion tensor, which characterize the polymer stretch, are plotted 
in Fig. 4b. We observe that polymer stretching is much more 
pronounced for planar and equibiaxial elongation compared 
to in planar shear flow. We have further validated the pro- 
posed algorithm for the rate (7) and chain length (N) depen- 
dence of the shear viscosity (see Fig. 5, which offers a quan- 
titative comparison with available NEMD data from [14] for 
an identical system). Since our method does not require flow- 
adapted boundary conditions, we are able to include here the 
first simulation results on steady state equibiaxial elongation. 
All results for the sample application, including the many be- 
yond the scope of this article and therefore not reported here, 
reproduce available experimental findings for gyration ten- 
sor and viscosities (shear thinning, strain hardening only in 
simple elongation, alignment in shear weaker than in elon- 
gation at same flow invariants, scaling behavior, overshoots, 
cf. [2, 27, 28]). The results provided in this section clearly 
demonstrate that the proposed simple procedure outlined in 
Tab. I allows to (i) recover known results obtained via classi- 
cal approaches, (ii) study flow geometries not accessible using 
alternate approaches, (iii) calculate the friction matrix and La- 
grange multiplier, i.e., the irreversible part of the closed and 
low-dimensional time evolution equation (5) for the coarse- 
grained variable in a straightforward manner. 



V. CONCLUSIONS 

Using an alternating MC-MD iteration scheme, our ap- 
proach successfully bridges the time-scale gap between mi- 
croscopic and macroscopic scales by establishing the coarse- 
grained model within a nonequilibrium thermodynamics 
framework. Since only short MD simulations are needed, our 
method is very efficient (moreover, it is ideally suited for par- 
allelization) and particularly allows to deal with arbitrary flow 
gradients, since neither special boundary conditions nor other 
constraints are needed. To be specific, even from the view- 
point of material property determination, our method is more 
efficient than standard NEMD when 7T < (no/n)(r/r s ), 
where ?io denotes the number of strain units needed for the 
NEMD. With t/t s = 30, n = 500 used here, taking n = 10 
from [15] and also the time for the MC step into account (see 
[23]), our method is superior to NEMD for 7T < 0.5, and m 
orders of magnitude faster at a value 10™ times smaller than 
that dimensionless rate. 

The presented approach is very general, but the (a- 
posteriori validated) success of the coarse-graining procedure 
depends crucially on the proper choice of the slow variable. 
As mentioned above in our illustrating example, conforma- 
tion tensors as slow variables for polymer melts are clearly re- 
stricted to the unentangled regime because interchain effects, 
entanglements or knots, hinder the relaxation of the confor- 
mation tensor for high molecular weight polymers [14, 16]. 
Some promising candidates for other soft matter systems are 
the tensorial order parameter for liquid crystals, the magne- 
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tization for magnetic liquids, and the path length of the en- 
tanglement network for entangled polymer melts [2, 14, 29]. 
The MC step is particularly challenging for dense polymeric 
systems, but efficient schemes exist for FENE as well as for 
atomistic models [26, 30]. 

For many complex fluids, Eq. (3) is known to serve as a suc- 
cessful starting point to derive closure relationships [14, 19]. 
Therefore, our method establishes the coarse-grained model 
all the way from equilibrium up to the validity of (3) and 
complements standard NEMD methods, which remain often 
well-suited for the less challenging regime of strong external 
forcing. 
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out of equilibrium like in a flow situation, several of the low- 
est Rouse modes are typically excited. In order to address 
this issue, we propose to include all Rouse modes in a single 
quantity such that the increasing relaxation times of the higher 
modes are reflected in a decreasing weight c p . Such a choice 
can be motivated by the fact that, in a stationary flow situation, 
the Lagrange multiplier is proportional to the product of relax- 
ation time and velocity gradient (see e.g. Eq. (8.52) in [13]). 
Thus, we use a single Lagrange multiplier A in the nonequi- 
librium ensemble Eq. (3) in order to excite all Rouse modes 
at the same time in a way that is consistent with Rouse theory. 
With c p = N/(irp) 2 decreasing for increasing mode number 
p as the corresponding relaxation times, x becomes the gyra- 
tion tensor, at least to a very good approximation. Comparing 
different choices for x (gyration tensor and the second rank 
tensor formed by either the first Rouse mode or the end-to- 
end vector), we found the gyration tensor to give the most 
accurate results. 



APPENDIX A: CHOICE OF SLOW VARIABLES 



APPENDIX B: REWEIGHTING SCHEME 



The proper choice of appropriate slow (collective) vari- 
ables is crucial not only for the method proposed here but 
for a broad class of nonequilibrium statistical mechanics ap- 
proaches based on projection operator techniques [13]. 

In the present case of unentangled polymer melts, there 
is ample evidence that single chain conformation tensors are 
promising candidates for slow variables [2, 13, 26]. There- 
fore, we assume the slow variable x can be decomposed into 
an average over single chain (symmetric second rank) confor- 
mation tensors, 



(Al) 



a=l 



The latter can always be expanded in a series of Rouse modes, 

AT-l JV-1 

xW^c p xW X W; xW = ^n pj Qf, (A2) 

p=l j=l 



where Xp°' is the p-th Rouse mode of chain a, Q T 

\/2/N sm(pnj /N) an element of the Rouse matrix and 
the connector vector of particles j + 1 and j of chain a [2]. 

As possible choice of the weights c p in Eq. (A2), we ini- 
tially implemented c p = 5 p i, i.e. only the first Rouse mode 
is included. This choice is reasonable since the first Rouse 
mode is the slowest and therefore a natural candidate for the 
slow variable x. However, the resulting model is restricted 
to very small deviations from equilibrium because there is no 
clear time scale separation to the higher modes which are ne- 
glected. In fact, the relaxation time of mode p in the Rouse 
model is t p = (/{8ksm 2 (pn/2N)] oc Ti/p 2 (£ and k are the 
bead friction coefficient and the spring constant in the Rouse 
model, respectively) and therefore the second mode relaxes 
only a factor four faster than the first one. Thus, when driven 



We here describe the reweighting scheme employed in the 
above described algorithm. The generalized canonical dis- 
tribution (3) is denoted by pa(z), in order to make explicit 
its dependence on the values of the Lagrange multipliers A. 
Corresponding averages of phase space functions are (A) a = 
Jdz A(z)pa(z). The analytical form of (3) allows to relate 
the distribution pa+sa (z) corresponding to different Lagrange 
multipliers A + SA to p\(z) by 



/Oa+5a(z) = Pa(z) 



e -5A:n(z) 



(Bl) 



Therefore, also the averages of phase space functions corre- 
sponding to different values of Lagrange multipliers are re- 
lated by 



(A) 



A+5A 



(B2) 



For small deviations i5A, the latter expression simplifies to 



A+5A 



(A) A - 6 A : (IL4) 7 
1-SA: (n) A 



(B3) 



Equation (B2) or (B3) for A = II and A = A4 are used in the 
reweighting scheme in order to calculate the corrected values 
for x and M, respectively, from recorded averages. In princi- 
ple, Eq. (B2) allows to recalculate averages of A for arbitrary 
5 A. In practice, however, due to the finite ensemble size, such 
estimates are accurate only if p\ and pa+sa have consider- 
able overlap. This is the case for SA small enough such that 
relevant states for averages at A + 5 A are sufficiently well 
sampled with p\. 
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APPENDIX C: STRESS TENSOR IN GENERALIZED 
CANONICAL ENSEMBLE 

Point of departure is the microscopic expression for the to- 
tal stress tensor [2] which can be inferred from the term of sec- 
ond order in the expansion of the configurational Helmholtz 
free energy with respect to the strain tensor [31] 



'a/3 




(Cl) 



The first, kinetic contribution can be well approximated by 
the ideal gas expression with p = nk^T. Deviations from 
this expression are minor in polymer melts and show up only 
at extremely high flow rates [14, 32]. 

We further assume (i) potential forces Fa a = —dH/drj a , 
and (ii) a generalized canonical ensemble p(z) = 
(1/Z*)exp[-/3H - N ch A: LT], cf. Eq. (3), where iV ch de- 
notes the number of polymer chains. 

Assumptions (i) and (ii) allow us to write 



dH 
1 d 



P 9r 3,P 

Inserting this into (Cl) gives 



— 5-A M „- p. 



(C2) 



'a/3 



/3V ^ 



dz : 



3,a 



Or 



'3,0 



-P(») 



an 



//// 



3,a 



Or 



V(3 

r 3, a 



(C3) 



For the special case of conformation tensor models, LT can 
be expressed as a bilinear form of the particle positions. Then, 
we obtain from Eq. (C3) the final expression 



er tot = -pi - 2n p k B Tx ■ A. 



(C4) 



Equation (C4) can independently be derived from nonequi- 
librium thermodynamics [13]. It should be noted that Eq. (C4) 
captures the polymer contribution to the stress tensor as 
the Lagrange multipliers A describe nonequilibrium polymer 
configurations. For short-chain polymer melts, a "simple fluid 
contribution" has to be added in order to account for the stress 
contribution to the total stress tensor that would be present in 
the absence of chain connectivity [21, 33]. 
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